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Abstract 



The structure of the thermal equilibrium state of a weakly interact- 
ing Bose gas is of current interest. We calculate the density matrix of 
that state in two ways. The most effective method, in terms of yield- 
ing a simple, explicit answer, is to construct a generating function 
within the traditional framework of quantum statistical mechanics. 
The alternative method, arguably more interesting, is to construct 
the thermal state as a vector state in an artificial system with twice 
as many degrees of freedom. It is well known that this construction 
has an actual physical realization in the quantum thermodynamics of 
black holes, where the added degrees of freedom correspond to the 
second sheet of the Kruskal manifold and the thermal vector state is 
a state of the Unruh or the Hartle-Hawking type. What is unusual 
about the present work is that the Bogolubov transformation used 
to construct the thermal state combines in a rather symmetrical way 
with Bogolubov's original transformation of the same form, used to 
implement the interaction of the nonideal gas in linear approximation. 
In addition to providing a density matrix, the method makes it possi- 
ble to calculate efficiently certain expectation values directly in terms 
of the thermal vector state of the doubled system. 



1 Introduction 

The temperature of black holes, predicted and investigated by Jacob Beken- 
stein on the basis of classical thermodynamic and statistical reasoning [I], H , 
is intimately related to the mixing of normal modes in quantum field the- 
ory. In the work of Hawking || (see also @]) the disappearance of half the 
degrees of freedom of a field into the black hole was shown to be respon- 
sible for the entropy and temperature of the remaining degrees of freedom 
as an astrophysical black hole evaporates. Unruh [|J constructed radiating 
and (implicitly — cf. || [7]]) equilibrium states of a quantum field on the 
maximal Kruskal extension ( "eternal" black hole) in terms of coherent com- 
binations of normal modes on the two sheets of the wormhole. The analytic 
continuation involved in the Unruh construction (which tells precisely how 
to combine positive-frequency modes on the physical sheet with negative- 
frequency modes on the unphysical sheet to create a state for the physical 
modes alone that is thermal at infinity) was soon recognized 0, ||, [UJ to be 
an instance of the generic analyticity property of thermal states in a complex 
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time coordinate [ITT], [T^|. Parker [O] showed that the well-known mixing of 



positive- and negative-frequency modes in cosmological models would give 
rise to a thermal spectrum of created particles under some (restricted but 
plausible) conditions. (In Parker's scenario the two members of a correlated 
pair are both physical, but they are effectively rendered decoherent by spatial 
separation.) 

Israel |14| and Sewell recognized Unruh's construction as a physical 
realization of a more abstract construction already known in quantum sta- 



tistical mechanics [15], [12], [16| . In the latter, a fictitious copy of the physical 
system is introduced and its modes are mixed with the physical modes to 
produce a thermal state of the latter. (We review the details in the next 
section, for the case of a free boson system.) In Unruh's scenario the second 
set of modes is not fictitious. At worst, it resides on the second sheet of the 
Kruskal wormhole; for a uniformly accelerating observer in flat space-time 
or an inertial observer in de Sitter space, it belongs to the part of ordinary 
space-time beyond the observer's horizon [|, [L7|, [7j . 

Black-hole temperature and cosmological particle production are effects 
in linear (noninteracting) quantum field theory in curved space-time that 
amount mathematically to linear mixing of creation and annihilation opera- 
tors in a way that preserves their canonical commutation (or anticommuta- 
tion) relations. Typically one set of operators defines a Fock-space structure 
so that the state under investigation is the corresponding vacuum, or no- 
quantum, state, while the other set of operators corresponds directly to phys- 
ical observations. This type of construction (which pervades |3|, || |5], [7], p~3j| ) 
is called a Bogolubov transformation. 

It is important to note that originally the Bogolubov transformation had 
nothing to do with either gravitational physics or temperature. N. N. Bogo- 
lubov [|T^, ^{J introduced it to find the ground state of a system of interacting 
massive bosons, such as atoms. More precisely, the original application was 
to ultracold helium, for which it turned out not to be quite adequate; but 
nevertheless, it has become a central tool in condensed-matter physics. In 
recent years, the necessity of a more detailed understanding of ultracold 
gases of weakly interacting bosons brought about many new applications of 
Bogolubov's formalism in physical contexts where the underlying approxima- 
tions are quite reliable. A recent article of Zagrebnov and Bru [^1J discusses 
these approximations, provides additional references to Bogolubov's pioneer- 
ing work, and reviews the current mathematical status of the theory. 

Our exposition of Bogolubov's theory follows that of Kocharovsky et al. 
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]. Under the assumption of a two-body interaction, the Hamiltonian is of 
the form 

h 2 k 2 

H = ~2M a ^ aVi + ^ t/ kik2,k 3 k4 a ki a k 2 a k :i ak4 • (1) 
k kj 

If the interaction U is sufficiently weak, most of the atoms will remain in 
the ground mode of the noninteracting theory ("Bose condensate"), and the 
excitations out of that mode can be treated to good approximation by a 
linearization of H. One introduces new operators 

A)= , a ° , Ac = /?oVc for MO. (2) 
'a Q a + 1 

The operators with k / satisfy the canonical relations [(3^, = Sy_y, 
and one can assume that terms in H containing more than two of them are 
negligible, so that H is approximately of the form 

H B = E + J2 VkPlfa + E(^WU + H.c), (3) 

k^O k^O 

with coefficients independent of the sign of k. (Conservation of momentum 
and of parity have been assumed here. At the next step we assume positiv- 
ity of the energy.) It is well known how to "diagonalize" such a quadratic 
Hamiltonian by a Bogolubov transformation: 

Ac = u k b k + v k &t k , K| 2 - |f k | 2 = 1; (4) 
H B = E + ^ b %- (5) 

k^O 

(All the undefined quantities in (|3J), (H), (^J) can be calculated from the U and 
M in (JJ); the details are not important here.) The dynamics of the system 
in the b variables is now trivial; for physical interpretation one returns to the 
(3 variables, for example, to see the structure of the physical ground state in 
terms of real atoms. 

In the wake of the recent experimental observations of Bose condensa- 



tion of atoms [23], interest in the type of problem originally investigated by 
Bogolubov has intensified. Methods have been invented to measure directly 
the (magnitudes of the) amplitudes «k and vy_ in ((D in a Bose-Einstein con- 



densate p3. In [E3] the (weakly) interacting Bose gas was studied at (low 



but) finite temperature; formulas were found for the characteristic function 
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(Fourier-transformed probability distribution) of the total number of atoms 
in a pair of modes, /?k/3 k + /3l k /3_k , and some related statistical quantities 
(the generating cumulants). In the interaction was treated as above, and 
the temperature was treated by constructing the statistical operator (density 
matrix) for each mode in the traditional, direct way: 

p k = (1 _ e -^) e ^/r. ( 6 ) 

(Note that it is the dressed (b) quanta that appear in (|), because the Hamil- 
tonian is (||), but it is the bare (J3) quanta whose statistics are being studied.) 



The formulas of [22[ are expressed in terms of a quantity 



and z(—Ak), where = t>k/wk is the parameter defining the Bogolubov 
transformation (H). 

We recall that when and t>k are real-valued, which one may assure 
by choosing the phases of the modes properly, it is customary to write the 
coefficients in (f|) as 

«k = cosh^k, fk = sinh^k! A^ = tanh^k . (8) 

(Henceforth we drop the subscript k whenever no ambiguity results.) Then 
9 is an additive quantity precisely analogous to the rapidity parameter of 
a Lorentz transformation; that is, the result of successively applying two 
Bogolubov transformations is the Bogolubov transformation corresponding 
to the sum of the respective 9s. Now (|7|) superficially resembles the addition 
formula for the hyperbolic tangent function, 

tanhx ± tanhy 

tanh(x ± y) = — — . (9) 

1 ± tanhxtanhy 



This suggests that the results of [g2| can be better understood, and perhaps 
more easily derived, in terms of the composition of the transformation of the 
original Bogolubov type (implementing the interaction) with a Bogolubov 
transformation of the Unruh type to implement the thermalization. Unfor- 
tunately, there are three things wrong with this naive idea. First, e e//T is 
greater than unity, and hence cannot be a hyperbolic tangent. Fortunately, 
that problem is instantly solved by rewriting (|7|) as 

1 Ae t/T - 1 A- e~ e/T 

(10) 



z(A) e"/ T -A 1 - Ae-/ T 
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which really does have the form (|9|). Second [see fl2"0|)1, the correct Bogo- 
lubov parameter in the Unruh (more properly, the Araki-Woods-Takahashi- 
Umezawa) construction is e _<E//2T , not e~ t ^ T . And third, the two Bogolubov 
transformations are about "different axes" in a space of dimension greater 
than 2, so there is no reason to expect them to compose in such a simple 
way. 

' " ™ in 



The present authors rederived (and generalized) the results of |22 
a simpler way, without, however, elucidating any connection with the thermal 
Bogolubov construction. Here we present (in Sees. || and §) a direct assault 
on the problem by the thermal Bogolubov method. Inevitably, the treatment 
displays a certain "hidden symmetry" between the parameters 



6 = tanh" 1 A and = tanh" 1 e" e/2T . (11) 

Equally inevitably, this symmetry is largely destroyed in the final formula for 
the density matrix, since to get it one must trace over the fictitious modes of 
the thermal construction but does not trace over any of the physical modes 
involved in the Bogolubov dressing transformation (see Sec. |||). However, the 
symmetry is restored when one takes the expectation value of an observable 
that involves only one physical traveling-wave mode. In fact, the moments 
and cumulants associated with the number operators for one or two modes 
can be most easily calculated directly from the Araki-Takahashi-Unruh pure 
state vector, rather than from the density matrix (Sec. |5|). In Sec. to obtain 
a simpler form for the density-matrix elements (generalizing p5[) we employ 
a different method, based on generating functions. 

2 Main calculation 

The analysis proceeds in seven steps. 

2.1 Reduction to standing- wave modes 

Introduce (for each pair ±k) 

K= —?r< b ' = —jr< (12) 

so that 

6 k e ik x + 6_ k e _ik x = V2[b c cos (k-x) + ib s sin (k-x)] , 
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blh + 6l k 6_ k = b% + b\b 



(13) 



and these operators and their adjoints satisfy the canonical commutation 
relations. Henceforth we concentrate largely on the cosine modes and often 
omit the subscript c. 



2.2 Dressing by the interaction 

The Bogolubov transformation (^) and (|8]) is chosen to convert (|3]) into ([5]). 
The inverse of such a transformation is the same except for the sign of 9, A, 
and v. The transformation thereby induced on the cosine modes is 

(3 = cosh9b + smh9b j ; b = cosh 0/3 - sinh 9ft; 9 = tanh -1 A. (14) 

The transformation of the sine modes has the opposite sense: 

(3 S = cosh 9 b s - sinh 9b\; b s = cosh 9 (3 S + sinh 9 (3\. (15) 



(This strange sign is an artifact of our phase convention in (|12|); it could be 
avoided at the cost of unnecessary imaginary numbers elsewhere.) 

Recall that f3 annihilates bare quasiparticles, here interpreted as excited 
physical atoms, while b annihilates dressed quanta, so that the eigenstates of 
energy are eigenstates of the number operator b^b. If |0) is the ground state: 
6|0) = 0, and \n) is the state of n physical atoms: ft/3\n) = n\n), then 



(2i - IV 

\2j) (A = tanh#). (16) 



It is instructive to review the proof of flTBD : Write |0) as iV J2^=o c n \n) with 
Co = 1. Impose the condition that (cosh#/3 — sinh#/^)|0) = to obtain the 
recursion c n+ i = J^h tanh 9 c n _ x . It follows that c n = for odd n and c 2 j = 



A 3 \J (2j — l)!!/(2j)H (where nil = n{n — 2) • ■ ■). Impose the requirement that 
(0 1 0) = 1 and observe that 

1 (17) 



to conclude that N 2 = y 1 — tanh 9 = (cosh#) 1 . 
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2.3 Thermal dressing 

The thermal Bogolubov construction must be applied to the energy eigen- 
states, regardless of whether there is an interaction. So (working with one 
mode at a time) we start from the number eigenstates satisfying 

6|0) = and b ] b\n) = n\n) . (18) 

Now introduce a fictitious mode "on the other side of the world" with cor- 
responding operators b and 1j; we work in a doubled Fock space with basis 
vectors satisfying 

b\Q, 0) = , 6|0,0) = 0; b^b\n, m) — n\n, m) , b b\n, m) — m\n, m) . 

(19) 

The thermal Bogolubov transformation is [16|, || 
c = cosh 05 — sinh (plj , c = cosh06 — sinh06^; e~ e ^ 2T = tanh 0, (20) 

with inverse b = cosh c + sinh <f> c>, b = cosh c + sinh c*. 

Let |0, 0} be the state with no c and c quanta: c|0, 0} = = c|0, 0}. Then 

-| oo oo 

|0, 0} = — — tanh n \n, n) = (1 - e~ £/T ) 1/2 £ e - ne/2t \n, n). (21) 
cosh n=0 n=0 

The proof is parallel to that of (|1(J) and slightly simpler, because the alter- 
nating-factorial and consequent square- root functions do not arise. 

The statistical (density) operator for the whole system (of one real and 
one fictitious mode) is 

oo 

|0,0}{0,0| = (l-e" e/T ) e-( n+m)e/2T \n,n)(m,m\. (22) 

n,m=0 

The statistical operator for the physical mode alone is obtained by tracing 
over the states of the unphysical mode: 

oo 

p = {l-e-^ T )Y j e- ne/T \n){n\. (23) 

n=0 

By design, p is precisely [cf. (||)] the thermal ensemble at temperature T! 
This is the generic Araki-Takahashi-Unruh construction, in the setting of 
our particular problem. 
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The foregoing notation is appropriate for the cosine modes. For the sine 
modes, we choose to change the sign of in ( p0|) . This convention cor- 
responds to the natural, momentum-conserving Bogolubov transformation 
(|43| ) on the original traveling-wave modes (where 6k mixes with &l k > n °t 
6 k ), which is mandatory in an Unruh-type situation but arbitrary when the 
barred modes are completely fictitious. The sign cancels in the final reduced 
density matrix, (|23|). 



2.4 Diagonalization of the composite transformation 

We now combine (0) and ( pop to express the bare modes in terms of the 
doubly dressed modes: 



(P \ 

P 

\P ] J 



I up vp vq uq \ 
vp up uq vq 
vq uq up vp 













c 


) 





(24) 



(Here p = cosh 0, q = sinh 0, and j3 is related to b just as (3 is to b.) Introduce 



n = e + 



and 



G 



P + P 



H 



P-P 



y/2 ' y/2 

Then a short calculation shows that 



* = e 
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c + c 



5 



c — c 



(25) 



(26) 



(27) 



That is, ( p4|) decouples into two elementary Bogolubov transformations! 
Moreover, the idea that © should be playing a role is vindicated by (f2~5|). 



y/2 ' v/2 
G = cosh 7 + sinh fi 7 f , if = cosh ^ S + sinh $ 



2.5 Construction of the thermal state with interaction 

Note that the state annihilated by 7 and 5 is the same as |0,0}, the one 
annihilated by c and c; after tracing over the barred mode, it will give us the 
thermal state we want. But to interpret that state we need to express it in 
terms of the /3'/3 number observable. 
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The first step is to construct |0, 0} within the (G, H) Fock space. This is 
done by applying the mathematics of ([16]) to the two transformations (j27]) 
independently. We need to define yet another basis: 

G^G\J,L] — J\J,L], H ] H\J,L] = L\J,L]. (28) 

Then 



\2j,2l}. 
(29) 



|0,0} = 1 V tanh J fttanh^ (2j . ^-^L ~ 

1 1 Vcoshficosh^^ \| (2j)!!(2/)!! 

The basis of physical interest is defined by 

P^(3\n, m) — n\n, m), /3 /3\n,m) = m\n,m). (30) 

In view of fl26|), the connection between \J,L] and \n,m) is just a unitary 
transformation, introducing no further pair creation. (Nevertheless, it is the 
source of most of the combinatorial complexity of our result.) One has 

\2j,2l) = —^(G^-^(HY\0,0] 

2^ l y/(2j)\(2l)\^ f \ K J W 

x (ft) 2j - K+2l - x 0) K+x \Q,O) 
1 y, (-l)V(2j)!(2Q! 



" 2i+ l k! (2j - «)!A! (2/ - A)! 

x ^/(2j + 2/ - k - A)! (k + \)\\2j + 21 - « - A, k + A). (31) 
Introduce m = k + A and define 

min(m,2Z) (~ 1) A 

Q{m ' M = k-S^fl A!(2(-A)!(m-A)!(2.,- m + A)! ' (32) 

Then 



/(2?)' (2/)' 2 i +21 

\2j,2l] = VV ^ Y, ^^J + 2l-m)\m\Q(m,jJ)\2j + 2l-m,m). 



m=0 

(33) 
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Substituting fl3H| ) into (|2"5| ) and simplifying, one obtains (for a cosine 
mode) 



1 



|0,0} = —j==== ^2 \ n i Tn)y/n\ ml 
V cosh f2 cosh \1> n m=0 



n+m even 

P 

X 



V tanh p "' Q tanh' ( 2/ ) ! ( 2 P 2 ! q ( _ j ^ (34) 
^ 4?ll(p-l)l v y v ; 



where p = |(n + m). For a sine mode, a factor (— l) p should be inserted. 



2.6 Evaluation of Q 

First note that whenever A is outside the range specified in (0), the summand 
is because at least one of the denominator factors is at a pole of the gamma 
function. Therefore, one may extend the summation over — oo < A < oo, 



and it is not necessary to write the limits at all |26 
Now define 



Pfa, c, m) 



A 



m 



Then, on the one hand, 



P(a,c,m) =P^- a - m ' a - m \0), 



(35) 



(36) 



where P^'^ (x) is a Jacobi polynomial [[27] . (This special value of the Jacobi 
polynomial apparently cannot be reduced to anything simpler; the very sys- 
tematic software associated with [28| identifies it as a certain hypergeometric 
function (also in [[27[) but nothing less.) On the other hand, we have 

4~ p (2l)l (2p- 21)1 Q(m,p -1,1) = 2~ n P(2l 1 2p,m) (p=\{n + m)). (37) 

In this notation (^) becomes 

1 



|0,0} 



V cosh Q cosh \l/ 



E 



n,m=0 
n+m even 



X 



y tanh p ^ tanh' ^ — 

On 



n, m)Vnl ml 

P(2l,2p,m). (38) 



1=0 



2 n ll(p-l)l 
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2.7 The density matrix 



The density matrix p nn i = (n\p\n f ) for a cosine mode is obtained by tracing 
|0, 0}{0, 0| over the states of the fictitious partner mode. That is, in a linear 
combination of objects \n,m)(n' ,m'\ one must set ml — m and sum over m. 
Since n — m and n' — w! are constrained to be even, the result is unless 
n — n! is even. One gets 

Pnn 1 = if n — n' is odd, (39) 



1 y/n\ n'\ 

cosh Q cosh ^ 2 n+n ' 



E 



ml 



m=0 
n—m even 



X 



E E P( ' / ;, 2 r m) ^' r,V ;r ) tanh— ' a tanh^' *, 



n + m . n' + m , . . 

p = — - — , p = — - — , n — n even. (4(JJ 

(The limits on the I and V summations are superfluous, for the same reason 
explained earlier for A.) 

The density matrix for a sine mode is the same except for a factor 
( — l)( n ~ n ')/ 2 . The density matrix for the entire atomic system is the ten- 
sor product of all the density matrices for the individual modes. (Recall that 
the latter depend on |k| through Q and 



3 Symmetries (or not) 



The Jacobi number (|35|) has these symmetries: 

P(c — a, c, m) = (— l) m P(a, c, m) ; 
2 m P(a,c,m) = (-l) a 2 n P (a, c,n) if n + m = c. (41) 

The first of these (with c = 2p, a = 21) expresses the essential invariance of 
our formulas when a summation index I is changed to p — I, hence when Q 
interchanged with The second (with a even) expresses the symmetry of 
|0, 0} in the real and fictitious modes. 

It is time to contemplate the meaning of Q and \& in terms of the basic 
parameters of our problem, A and e/T. From the definitions we have 

A + e~ t/2T A - e~ e/2T 

tanhn= ^ tanh^= 1 _ Ae _ e/2r . (42) 
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These bear tantalizing resemblances to z( Z A ^ and —j^) as ex P resse d through 
(|T0|), but the factor of 2 with the temperature is ineradicable. 

Now observe the behavior of the density matrix (|40|) under elementary 
operations on the parameters: 

• Interchanging Q and \& (i.e., changing the sign of 0) leaves p invariant. 
Since this operation amounts to replacing e~ € ^ 2T by a negative number, 
it may appear unphysical. However, as previously remarked in connec- 
tion with the sine modes, it really represents just an arbitrariness in 
the definition of the fictitious modes. 

• Changing the signs of both Q and \1/ is equivalent to changing the signs 
of both 8 and (p. Its effect on p nn i is an overall factor (-l)("-»')/ 2 . This 
is precisely the distinction between cosine and sine modes. 

• Interchanging 9 and is equivalent to changing the sign of Its ef- 
fect on p is substantive: Each term in the summand is multiplied by 

Thus (of course) the final formulas of the theory are not all 
symmetric under interchange of an interaction parameter with a tem- 
perature parameter, despite the intriguing symmetries in the formalism. 
(We shall see later, however, that some formulas are symmetric.) 

The other 4 nontrivial elements of the group generated by these operations 
add no additional insight. 



4 A pair of modes 

At the cost of dealing with twice as many indices, one can work directly 
with the original traveling- wave modes, skipping step ([121) . Some of the 
intermediate results are useful, so we summarize the calculation here. 
Define a momentum-conserving thermal Bogolubov transformation 

bk = PkC k + <?kcl k , & k = £>kC k + gkcL k (43) 

equivalent to (|^) and its b s counterpart. When ( f4~3"D is combined with (||) 
and its barred counterpart, one obtains two 4x4 systems of precisely the 
form (|24"D , except that all the creation operators have the opposite sign of k 
from the annihilation operators they mix with. So one can define operators 



as in (|26|), with subscripts ±, and get 



G± = cosh Q 7± + sinh Q 7^ , H± = cosh ^ 5± + sinh * 5^ . (44) 
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By the same methods as before, one finds that the analog of (P9"|) is 



'°' ' 0,0 ' = cosh » cosh* E 'anh^tanh'*!,,,] <8 |/,/], (45) 

3,1=0 

where |0, 0, 0, 0} is annihilated by y± and 5± , and 

&+G+M=M5\, HlH + \lJ] = l\lJ], 
GtG-\j,3\=Mj\, HlH_\l,l] = (46) 

Note that in fl45|) the occupation numbers may be either even or odd, but 
the number of G + quanta is always the same as the number of G_ quanta, 
and similarly for H. The analog of flUf ) for standing waves is the tensor 
product of (|29|) for the cosine mode with its partner for the sine mode; in 
that case, all occupation numbers must be even, but there is no constraint 
relating those for c quanta to those for s quanta. 

Let p n h,n'n' be the matrix element of the statistical operator between the 
state (nn\ with n (physical) (3 + quanta and n /3_ quanta and another such 
state \n'n'). A long calculation parallel to that leading to (HO) yields 



Pnnn'n' =0 if U ~ U ^ Tl' - ft' , (47) 

and otherwise 



1 ^2, ^ X, \fn\ h\ n'\ n'\ ml ml 

Pnfinlfl ' = cosh 2 ft cosh 2 ^ ^ ^ 2f+f'- 2m - 2 - II VI (p-l)l (p> - V)\ 
x P(l, p, m)P(l, p, m)P(l', p', m)P(l', p', m) 
x tanh p+p '-^ r n tanh l+l ' (48) 
m = n + m — n, p = n + m, p' = n' + m. 



A corollary of the momentum conservation constraint (f47|) is that n + n 
has the same parity as n' + n' in any nonvanishing matrix element; this 
is necessary for consistency with ([39]) and with the general principle that 
Bogolubov transformations "create" quanta only in pairs. 

5 Number observables 

Henceforth it is convenient to write formulas in terms of 

A = tanh6i and B = tanh</> = e" e/2T . (49) 
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Along with we have (noting that \AB\ < 1) 
cosh n = — = , sinh = 



1-A 2 )(1-S 2 ) ~ A 2 )(l -B 2 ) 

1 - AB A-B . . 

cosh * = . , smh * = . . (50) 

;i - A 2 )(l-B 2 ) ^ {1 - A 2 ){1 - B 2 ) 

Our formalism provides two ways to calculate the expectation value of 
an operator. One may use the density matrix, ([40"D or (|48|) . Alternatively, 
one can express the operator in terms of the operators G and H and their 



adjoints and take its matrix element in the pure state vector (29) or (45). 
For some observables the second method is much easier and also displays the 
6 <-> (or A <-> B) symmetry to the maximal extent. 

Consider a fixed mode pair ±k. Recall that all our creation and annihila- 
tion operators carry subscripts c, s, +, or — ; we routinely omit not only any 
reference to the vector k but also the subscript (especially c or +), when no 
ambiguity results. When the sine mode is involved, the correct state vector 
is the product of (^9|) with its sine partner, whose formula is the same as 
( |29| ) except for a factor (— in the summand. For any mode type there 
is a number operator 

n = fip = -—(G^ + H^) — (G + H) = -(G^G + H^H + G^H + H^G). (51) 
V 2 V 2 2 

5.1 The mean number 

First consider (n + ) or (n_). From (f4"5|) we see that only G^G and H^H 
contribute to (n) = {0, 0, 0, 0|n|0, 0, 0, 0}, because the other two terms in 
([51] ) destroy the equality of the + and — occupation numbers. Using the 
first of 

g^(T^r g^w^ 1 -*^ (52) 



along with (|50|), one calculates 

: ";i - a 2 )(i - B 2 ) 



I \ A2+B2 

\ n ) = n — i^r, — ™ • ( 53 ) 
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After chasing through several layers of definitions, one sees that (53) 



ought to coincide with (each term of) equation (72) of p2| , and it does. 
Incidentally, (72) is one of very few formulas in [|22[ that are (even implicitly) 
symmetric in A and B. But it is now clear why it, and any expectation value 
that concerns only one traveling-wave mode, must have that symmetry: The 
thermal Bogolubov transformation (^) connecting that mode to one of the 
fictitious barred modes has the same algebraic form as the transformation 
(|43| ) connecting that mode to another of the physical modes. When the 
expectation value is calculated, all other modes are effectively traced over, 
so there is nothing to distinguish the fictitious partner mode, with mixing 
constant B, from the physical one, with mixing constant A. 

This argument does not apply to standing waves, because the transforma- 
tion (|i4T) -(|i~5|) does not have the same form as (f43[) ; in fact, we shall soon see 
that the conclusion does not hold for such modes. Nevertheless, a calculation 
based on (|29|) shows that 

(n c > = (n s > = (n>. (54) 

In this case the second of formulas (0) is used, and again the G^H and H^G 
terms in (|5l| ) do not contribute, but this time for a different reason: they 



produce occupation numbers of odd parity, which do not occur in {0, 0| 

5.2 The second moment 

The operator n 2 is 



,r = \{O r + O p + O i ), (55) 



r 2 

where the terms 

R = (G^G) 2 + {H^Hf + AG ] GH ] H + G ] G + H^H (56) 
are relevant, the terms 

P = (G ] fH 2 + (H^fG 2 (57) 

are partially relevant, and the terms 

1 = 2[{G^) 2 GH + {H^fHG + G ] G 2 H ] + H^H 2 G^ + &H + H ] G] (58) 

are irrelevant because they contribute nothing to (n 2 ). The partially relevant 
terms conserve parity but not momentum, so they contribute in the context 
of (BSD but not (M). 
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For (nV) or (n 2 _) a calculation like that of (n) leads to a formula that can 
be abbreviated as 

(n 2 ) = 2(n) 2 + (n). (59) 

Since ((n — (n)) 2 ) = (n 2 ) — (n) 2 , it follows that the second moment of the 
number distribution is 

((n-(n)) 2 ) = (n) 2 + (n). (60) 

Again, these quantities are symmetric in A and B. 
A parallel calculation yields 

<„■> = <« 2 > = (i + fl 2 )[^ + a 2 + ^(i + fl 2 )] (61) 

Here the symmetry between A and B is destroyed by the contribution of the 
partially relevant terms. 

5.3 The second cumulant 



Kocharovsky et al. |E| do not give a formula for (n 2 ); the closest point of 
comparison is the generating cumulant k 2 . One case of their formula (71) is, 
after reexpression in terms of A and B rather than z(A) and z(—A), 

A 2 B i + A ± + 4 A 2 B 2 + B * + A 2 

K2 ~ (1-^)2(1-^2)2 ( 62 ) 

= <^> + A2(B2 + 1)2 



(1 - A 2 ) 2 {l-B 2 ) 2 ' 

Note that the second term is not symmetric in A and B. 

This cumulant refers to a pair of modes, ±k. Its definition is 

£ 2 = \\{(n + + n_-2{n)) 2 - 2(n)] 
= (n 2 ) - 2(n} 2 - (n) + (n + n_) 
= (n+n_), (63) 

where (B^) has been used at the last step (and (n±) = (n) from the begin- 
ning). Thus the asymmetry of in the two parameters arises from "in- 
terference" between two modes. The quantity is readily calculated by our 
method: 

n + n_ = \(0 N + T + ), (64) 
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where 



N = G ] GG ] _G_ + H ] HH ] _H_ + G ] GH ] _H_ + H^HGIG_ , 

T = G ] G ] _HH_ + GG_H ] H ] _ + G ] G_H ] _H + G ] _GH ] H_ , (65) 

and Oo comprises 8 terms with either three G operators and one H or vice 
versa. It is easy to see that Oo and the last two terms of Ot do not contribute 
to the expectation value of n + n_ in |0, 0,0,0}. The calculation of \{On) is 
very similar to that of (n 2 ) and results in the symmetric expression 



A 4 + 4A 2 B 2 + B 4 + \{A 2 + B 2 )(l + A 2 B 2 ) 
(1 - A 2 ) 2 (l -B 2 f 

For | (Ot) one finally gets an asymmetric expression, 

\{A 2 -B 2 ){1- A 2 B 2 ) 
(1 - A 2 ) 2 (l -B 2 ) 2 ' 



(66) 



(67) 



Adding (pq ) and (67) and rearranging, one reproduces 



5.4 The quadratic moments for standing waves 

The analog of (n + n_) for the standing modes is (n c n s ), the expectation value 
being taken in |0, 0} for the cosine and sine modes independently. Thus 



(n c n s ) = (n c )(n s ) = (nf 



A 2 + B 2 



1 2 



;i - a 2 )(i -5 2 ) 



Clearly one must have 

((n c + n s ) 2 } = ((n + + n_) 2 }, 

whence by (|^) and (|68D one has 

\[«) + (n 2 )} = -(n c n s ) + (n 2 ) + (n + n_) 
= (n) 2 + (n) + (n + n_), 



(68) 



(69) 



(70) 



which simplifies to the quantity ( pi]) upon substitution from (|5^) and 
fl63|). Given that (n 2 ) = (n 2 ), and that (n + nJ) is known, this argument is 
the easiest way to derive (|6l|) . 
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6 The generating function 

6.1 The single- mode density matrix revisited 



The matrix elements p nn > = (n\p\n') of ( p9|) and (^UJ) have been obtained 
by regarding the statistical operator p of (^3|) as referring to one mode of a 
two-mode system that is in the pure state |0, 0}{0, 0|. A different, and some- 
what more direct, method of calculating p nn , first constructs the generating 
function 

,(*,,)= t (71) 
n ,n'=o yn\n\ 

and then expands it in powers of x and y. The bra and ket states that show 
up here as the summations over n and n' are recongnized as the well known 
coherent states (of Bargmann type), the eigenbras of ft and eigenkets of /?, 

oo T n/„| oo I A n' 

£^(/?t_ x) = , (/?-i/)E^7=r = o. (72) 



=0 Vn! n , =0 vn'l 

These eigenvector equations are exploited upon writing p as a normally 
ordered function of ft and /3; that is, all fts stand to the left of all /3s. We 
begin with 



1 _ e -e/T )e -(e/T)6t6 = ^ _ g-e/T) eX p(-(l - eT^; ft) 
e |(^ 2 -^ (l _ e - e /T } exp ^_ (1 _ e -e/T )/3 t. pyW-fP)B j (73) 



where 

oo y-n\rn 
n=0 /t- 

denotes the basic ordered exponential function, for which the much used 
identity 

(1 _ A) 6t6 = e" A6t;b (75) 

is a familiar application. 

The ft; /5-ordered version of p, 
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with 

p -e/T 

A = 1 



(coshtf) 2 -e- 2e / T (sinh^) 2 ' 
0c/rrs sinh^ cosh^ 



(cosh#) 2 -e- 2e / T (sinh#) 



2 



, 1 _ e -e/T 

J\ 2 - M = , (77) 

/ (cosh^) 2 -e- 2e / T (sinh^) 2 



can now be found by a variety of methods. It suffices to verify that it responds 
correctly to infinitesimal changes of 9. 

Now, when using ( |76| ) in ([71]), the identities ( [72] ) permit the replacements 
(3* —>■ x, (3 — > y, and then the summation is elementary, 

± (78) 

n ,„'=o Vn\n'\ 

We thus arrive at the explicit form of the generating function (ffl|), 



g(x, y) = ^A 2 - n*n e ^ 2 e (1 - A )^e^ V . (79) 

Such a two-dimensional Gaussian is a typical generating function for lin- 
ear systems, since all time-transformation functions are of this form if the 
Heisenberg equations of motion are linear. A recent example is the paramet- 
ric oscillator investigated by Rashid and Mahmood p9| , who combine the 
Maclaurin series of the three exponential functions in (|79|) and get an answer 
where one summation is still to be done, somewhat like the situation in (40) 
above. 

Here we wish to put a different approach on record. It exploits the well 
known fundamental generating function for Bessel functions of integer order 
(see 9.1.41 and 9.1.5 in |0Ofl), 



e Ht-i/t) = t k J k {z)= >T t k (-l)^ kl - k) (z), (80) 

k=— oo k=— oo 

and one of the generating functions for Gegenbauer polynomials of index 
k + i (see 22.9.5 in |3C 



e-(z/2)- k J k (z) = p k] ^l iy (z 2 + u^C^ (u/V* + tf) , (81) 
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valid for k > 0. Jointly they amount to 

\z(t-l/t)+U _ fo-|H( ; ^(|t|+t)/ y/^gdfcl-fc) 

" k hM {Zt) ( ^ (1*1)1(21*1 + /)! 

x (z 2 + u 2 fc^ ] (u/V^+rf) . (82) 

The left-hand side thereof turns into the product of exponentials in ([79]) if 
one puts 

zt = fix 2 , —z/t = fi*y 2 , u — (1 — X)xy , (83) 

so that 

(z 2 + u 2 f 2 = v/(l-A) 2 -^ xy, \ = 1 ~ A (84) 

V^ + m 2 J (I - A) 2 - 



Note that the latter, which is the argument of the Gegenbauer polynomial in 
([82]), does not depend on x or y. All dependence on x and y is in the form 
of explicit powers. 

Upon putting things together, we arrive at 



e ±tix 2 +(l-\)xy+±ti*y 2 



9(x,y) 

k=—oo 1=0 

(2\k\)\ 

X 



j (V (1 - A) 2 - t**v x v) 



(|*|)! (2|A;| + 

x ( ; I , (85) 

! '(1 - A) 2 - n*n 

where n = \ k | + k + I is the total power of x and n' = \k\ — k + 1 is that of y, 
their difference n — n' = 2k and sum n + n' = 2\k\ + 21 being even. We write 

Tl — TV 

k = — - — , I = n < = min{n, n'} , 2\k\ + / = n> = max{n, n'} (86) 

to present (|55]) as a sum over n and n'. This enables us to identify the matrix 
element (n\p\n') = p nn > , with the outcome 

i 

(n> — n<)! (n<\ \ 2 



Pnn> = jt? ttt — r V A2 ~ 

(±(n> -n<))! \n>!/ v 
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x (/V2)2 ( "- n <V/2)s (n '~ n<) ((1 - A) 2 - fi*fiy n< 

x ci\ {n >- n<+1)) ( , X " A ) (87) 
"l-\) 2 -p*p) 

iin — n' is even, and p nn > = if n — n' is odd. In addition to the expressions 
of (ff7|) we also meet here 

2 e- 2 ^(coshg) 2 -(sinhg) 2 
(1_A) (cosh^- e ^/T (sinhe )2- ( 88 ) 

In fact, the result fl8~T|) is a bit too general for our purposes because we 
have the case of \x = p*, which simplifies matters somewhat. Expressed in 
terms of A = tanh^ and B = e _e//2T we have, for n — n' even, 



i i 

12 \ 2 



(n>-n<)! (n K \y ( I - A 1 , (] _ ^ 



(|(n>-n<))! U>! / I 1 - A 2 5 4 



x 



'1 A(1-B 4 ) \^ (W> " W<) f B 4 -A 2 \z n< 
2 1 - A 2 5 4 J I 1 - A 2 5 4 J 



x G (j'»>-"< +»> , t 1 -- 4 ^ 2 | . (89) 

[1 - A 2 B i )(B 4 - A 2 ) 



6.2 Diagonal terms 



The Gegenbauer polynomials of index | are the Legendre polynomials, 
C n 2 = P n , so that the diagonal matrix elements are given by 

1-A 



(n\p\n) = p nn = J\ 2 - p*p (1 - A) -p*p) 2 P 



7(1 -A) 2 -pp. 



(90) 

Upon recalling the familiar generating function for Legendre polynomials, 

OO 1 

y t n p n (x) = -. , (9i) 

^ 1 ; VI - 2tx + 1 2 ' 1 ; 

we obtain the corresponding generating function for p nn , 

A 2 - pfji 



n=0 



1 - 2(1 - X)e iu + [(1 - A) 2 - p*p]e 2iu 
22 



(92) 



which one recognizes to be identical with (26) in [33] when the differences in 
notation are taken into account. 



6.3 Comparison 

The dependence of (|89|) on A and B, although complicated, is nicely consis- 
tent with the observations in Sec. [| Only even powers of B appear; this is the 



symmetry of (40) under the interchange of VL and \F Similarly, the behavior 
under change of sign of Vt and \1/ is implemented by the factor AS n> ~ n< '' 2 . 

To establish direct contact between ( |89| ) and ( [40"D we would need to ex- 
pand 

1 - A 2 B 2 

p nn > cosh Q cosh ^ = p nn , — — — — (93) 

(1 - A ){1 - B ) 

in powers of 

, ^ A + B , , T A - B , 

tanh Q, = and tanh * = , (94) 

1 + AB 1 — AB K ' 

which looks like an unfairly difficult homework exercise. We have evaluated 
([89]) and (|40|) numerically for a variety of values of the parameters, always 
finding excellent agreement. 

Thus the method of the generating function has provided a formula, fl89"D , 
that is simpler than the one provided by the method of the thermal Bogo- 
lubov transformation, fl40|). (Numerically, the latter involves more computa- 
tion and requires truncation of an infinite series.) Nevertheless, the thermal 
Bogolubov method has given us an interesting and elegant way of looking 
at the problem, which illuminates its symmetries and may lead to further 
insights in the future. Furthermore, for observables of the sort studied in 
Sec. [|the calculations based on the thermal pure state vectors, (|29|) or (f45|), 
appear to be competitive with the conventional methods. 
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